/* 2010: Copyright by Massimo Bernaschi, Giorgio Parisi, Leonardo Parisi */
#include <time.h>

#include "timing.h"

#include "dictionary.h"
#include "iniparser.h"
#include "util.h"
#include "hsgcpu.h"


//-------------------------------------------------------------------------------------------------------
// Usage
//-------------------------------------------------------------------------------------------------------
void Usage(char *program) {
  fprintf(stderr,"Usage: %s -i inputfile [-r spin_file] [-o (over-relax only)] [-b (heat bath only)]\n",program);
}


//-------------------------------------------------------------------------------------------------------
// Init
//-------------------------------------------------------------------------------------------------------
void Init(void){

  if(sizeof(REAL)==sizeof(float)) {
    zero=0.0f;
    one=1.0f;
    two=2.0f;
    three=3.0f;
    zerofive=0.5f;
    zerotwentyfive=0.25f;
  } else if(sizeof(REAL)==sizeof(double)) {
    zero=0.0;
    one=1.0;
    two=2.0;
    three=3.0;
    zerofive=0.5;
    zerotwentyfive=0.25;
  } else {
    writelog(TRUE,APPLICATION_RC,"Invalid size of REAL: %d",sizeof(REAL));
  }

  zeroEne = 0.0; //MESSA SOLO PER IL CALCOLO DELL'ENERGIA IN DOUBLE COME PROVA

  Ra = (REAL*) Malloc(ms_s);
  Rb = (REAL*) Malloc(ms_s);
  Rc = (REAL*) Malloc(ms_s);
  Ba = (REAL*) Malloc(ms_s);
  Bb = (REAL*) Malloc(ms_s);
  Bc = (REAL*) Malloc(ms_s);

  Jpx = (REAL*) Malloc(ms_j);
  Jpy = (REAL*) Malloc(ms_j);
  Jpz = (REAL*) Malloc(ms_j);
  Jmx = (REAL*) Malloc(ms_j);
  Jmy = (REAL*) Malloc(ms_j);
  Jmz = (REAL*) Malloc(ms_j);

  Energy = (REAL*) Malloc(ms_s);

  sdbL = findlog2(side);
  sdbLs = findlog2(side*side);

  sdbLm1 = sdbL-1;  //log_2(L) -1 used to divide by L/2
  sdbLsm1 = sdbLs-1;  //log_2(L*L) -1 used to divide by (L*L)/2
  maskL = (side>>1)-1;
  maskL2 = ((side*side)>>1)-1;
  maskL3 = (V(side)>>1)-1;
  Ldb2 = side>>1;
  Lsdb2 = (side*side)>>1;

  nspin=V(side)/2;

  randNum=V(side);

  aRand = (REAL*) Malloc(randNum*sizeof(REAL));
  bRand = (REAL*) Malloc(randNum*sizeof(REAL));
  cRand = (REAL*) Malloc(randNum*sizeof(REAL));
  Rand = (REAL*) Malloc(randNum*sizeof(REAL));

  srand48(seed);

}


//-------------------------------------------------------------------------------------------------------
// findlog2
//-------------------------------------------------------------------------------------------------------
int findlog2(int n){

  int i=0;

  if(n<1 || n&1) { writelog(TRUE,FINDLOG_RC,"Invalid value %d. Must be > 0 and even\n",n); }

  while(n>0) {
    n>>=1;
    if(n>0) i++;
  }

  return i;
}


//-------------------------------------------------------------------------------------------------------
// nextPow2
//-------------------------------------------------------------------------------------------------------
unsigned int nextPow2( unsigned int x ) {
  --x;
  x |= x >> 1;
  x |= x >> 2;
  x |= x >> 4;
  x |= x >> 8;
  x |= x >> 16;
  return ++x;
}


//-------------------------------------------------------------------------------------------------------
// FreeHost
//-------------------------------------------------------------------------------------------------------
void FreeHost(void){

  Free(Ra);
  Free(Rb);
  Free(Rc);
  Free(Ba);
  Free(Bb);
  Free(Bc);

  Free(Jpx);
  Free(Jpy);
  Free(Jpz);
  Free(Jmx);
  Free(Jmy);
  Free(Jmz);

  Free(Energy);

  Free(aRand);
  Free(bRand);
  Free(cRand);
  Free(Rand);

}


//-------------------------------------------------------------------------------------------------------
// Dump_Spin
//-------------------------------------------------------------------------------------------------------
void Dump_Spin(int count){

  FILE* tempfile;
  char namefile[MAXSTRLEN];

  memset(namefile,'\0',MAXSTRLEN);

  if(count>=0) snprintf(namefile,sizeof(namefile),"spin_dump_%d",count);
  else snprintf(namefile,sizeof(namefile),"spin_dump_start");

  tempfile = Fopen(namefile,"w");

  int i; for(i = 0; i<ns_j; i++){
    if(sizeof(REAL)==sizeof(float)) {
      fprintf(tempfile,"%f %f %f\n%f %f %f\n",
              Ra[i],Rb[i],Rc[i],Ba[i],Bb[i],Bc[i]);
    } else {
      fprintf(tempfile,"%lf %lf %lf\n%lf %lf %lf\n",
              Ra[i],Rb[i],Rc[i],Ba[i],Bb[i],Bc[i]);
    }
  }

  fflush(tempfile);

  fclose(tempfile);

}


//-------------------------------------------------------------------------------------------------------
// Initialize_Spin
//-------------------------------------------------------------------------------------------------------
void Initialize_Spin(REAL* sa_a, REAL* sa_b, REAL* sa_c){

  REAL a, b, c;
  REAL invnorma, norma;

  int i; for(i = 0; i<nspin; i++){

    norma=two;

    while (norma>one){

      a=(RANDZEROONE)*two-one;
      b=(RANDZEROONE)*two-one;
      c=(RANDZEROONE)*two-one;

      norma=a*a+b*b+c*c;

    }

    invnorma=one/SQRT(norma);

    sa_a[i]=a*invnorma;
    sa_b[i]=b*invnorma;
    sa_c[i]=c*invnorma;

  }

}


//-------------------------------------------------------------------------------------------------------
// Read_Spin
//-------------------------------------------------------------------------------------------------------
void Read_Spin(char *inputfile, REAL* sr_a, REAL* sr_b, REAL* sr_c, REAL* sb_a, REAL* sb_b, REAL* sb_c) {

  REAL a, b, c;
  int j=0, k=0;
  FILE *fpi;
  fpi=Fopen(inputfile,"r");
  fprintf(stderr,"Reading spin from file %s\n",inputfile);

  int i; for(i = 0; i<V(side); i++){

    if(sizeof(REAL)==sizeof(float)) fscanf(fpi,"%f %f %f",&a,&b,&c);
    else fscanf(fpi,"%lf %lf %lf",&a,&b,&c);

    if(i&1) {
      sb_a[j]=a;
      sb_b[j]=b;
      sb_c[j]=c;
      j++;
    } else {
      sr_a[k]=a;
      sr_b[k]=b;
      sr_c[k]=c;
      k++;
    }
  }
  fclose(fpi);

}


//-------------------------------------------------------------------------------------------------------
// Initialize_J
//-------------------------------------------------------------------------------------------------------
void Initialize_J(void){

  int j;

  switch (init_type_j){

  case 0:{
    for (j=0;j<ns_j;j++) {
      Jpx[j] = ((RANDZEROONE)>zerofive) ? -1 : 1;
      Jpy[j] = ((RANDZEROONE)>zerofive) ? -1 : 1;
      Jpz[j] = ((RANDZEROONE)>zerofive) ? -1 : 1;
      Jmx[j] = ((RANDZEROONE)>zerofive) ? -1 : 1;
      Jmy[j] = ((RANDZEROONE)>zerofive) ? -1 : 1;
      Jmz[j] = ((RANDZEROONE)>zerofive) ? -1 : 1;
    }
    break;
  }

  case 1:{
    for (j=0;j<ns_j;j++) {
      Jpx[j] = Jpy[j] = Jpz[j] =
        Jmx[j] = Jmy[j] = Jmz[j] = one;
    }
    break;
  }

  case -1:{
    for (j=0;j<ns_j;j++) {
      Jpx[j] = Jpy[j] = Jpz[j] =
        Jmx[j] = Jmy[j] = Jmz[j] = -one;
    }
    break;
  }

  case 2:{ //Caso Gaussiano
    for (j=0;j<ns_j;j++) {
      Jpx[j] = rand_gauss();
      Jpy[j] = rand_gauss();
      Jpz[j] = rand_gauss();
      Jmx[j] = rand_gauss();
      Jmy[j] = rand_gauss();
      Jmz[j] = rand_gauss();
    }
    break;
  }

  default: {
    writelog(TRUE,APPLICATION_RC,"Invalid initialization for J: %d", init_type_j);
    break;
  }
  }

}


//-------------------------------------------------------------------------------------------------------
// rand_gauss
//-------------------------------------------------------------------------------------------------------
REAL rand_gauss(void){

  REAL r,cosr;
  REAL a=one-RANDZEROONE;
  if(a<EPSILON) { a=EPSILON; }

  r=SQRT(-2.*LOG(a));

  cosr=COS((2*M_PI)*(RANDZEROONE));

  return r*cosr;

}


//-------------------------------------------------------------------------------------------------------
// Parameters_Reading
//-------------------------------------------------------------------------------------------------------
void Parameters_Reading(char *inputfile){

  dictionary *ini;
  char key[MAXSTRLEN];

  if(inputfile==NULL) { writelog(TRUE,APPLICATION_RC,"No inputfile"); }

  ini = iniparser_load(inputfile);

  if(ini==NULL) { writelog(TRUE,APPLICATION_RC,"Cannot parse file: %s\n", inputfile); }

  READINTFI(side,"Side");
  READINTFI(init_type_j,"Initialization type for J");
  READINTFI(nsteps,"Number of steps");
  READINTFI(overrelaxiter,"Number of overrelaxation iterations");
  READINTFI(outfreq,"Output frequency");
  READINTFI(dumpfreq,"Dump configuration frequency");
  READREALFI(h,"Magnetic field");
  READREALFI(Beta,"Beta (1/T)");
  READINTFI(seed,"Random seed");

}


//-------------------------------------------------------------------------------------------------------
// Calculate_Sizes
//-------------------------------------------------------------------------------------------------------
void Calculate_Sizes(void){

  ns_j = (V(side)/2);

  ms_s = ns_j*sizeof(REAL);

  ms_j = ns_j*sizeof(REAL);

}


//-------------------------------------------------------------------------------------------------------
// Calculate_Sizes
//-------------------------------------------------------------------------------------------------------
void ComputeRandom(REAL* rand_ptr, int mode){


  switch (mode) {

  case ZEROONE:{
    int i; for(i=0; i<randNum; i++) rand_ptr[i]=RANDZEROONE;
    break;
  }

  case ONEONE:{
    int i; for(i=0; i<randNum; i++) rand_ptr[i]=(two*(RANDZEROONE))-one;
    break;
  }

  case ZEROFIVE:{
    int i; for(i=0; i<randNum; i++) rand_ptr[i]=RANDZEROONE-zerofive;
    break;
  }

  default: {
    memset(rand_ptr,0,randNum*sizeof(REAL));
    break;
  }
  }


}


//-------------------------------------------------------------------------------------------------------
// main
//-------------------------------------------------------------------------------------------------------
int main(int argc, char **argv){

  char *inputfile = NULL;
  char *po;
  int istep;
  //int nblocks_red, nthreads_red;
  double ttime=0.;
  REAL normene;
  int j;
  char *readspin=NULL;


  TIMER_DEF;

  int i; for(i = 1; i < argc; i++) {
    po = argv[i];
    if (*po++ == '-') {
      switch (*po++) {
      case '?':
        SKIPBLANK
          Usage(argv[0]);
        exit(OK);
        break;
      case 'v':
        SKIPBLANK
          verbose=TRUE;
        break;
      case 'i':
        SKIPBLANK
          inputfile=Strdup(po);
        break;
      case 'o':
        SKIPBLANK
          onlyOverrelaxation=TRUE;
        break;
      case 'r':
        SKIPBLANK
          readspin=Strdup(po);
        break;
      case 'b':
        SKIPBLANK
          onlyHeatBath=TRUE;
        break;

      default:
        Usage(argv[0]);
        exit(OK);
        break;
      }
    }
  }


  Parameters_Reading(inputfile);

  if (side<=0 || !pow_of_two(side)) {
    writelog(TRUE,APPLICATION_RC,"Invalid side: must be > 0 and power of two\n");
  }

  Calculate_Sizes();

  Init();

  if(readspin==NULL) {
    Initialize_Spin(Ra,Rb,Rc);
    Initialize_Spin(Ba,Bb,Bc);
    Dump_Spin(-1);
  } else {
    Read_Spin(readspin,Ra,Rb,Rc, Ba,Bb,Bc);
  }

  Initialize_J();

  normene=one/(three*V(side));

  Energy[0]=zero;

  ComputeHEnergy(Ra, Rb, Rc, Ba, Bb, Bc,Jpx, Jpy, Jpz, Jmx, Jmy, Jmz, Energy,TRUE);

  printf("Initial Energy per spin on host: %f\n", normene*Energy[0]);

  for(istep=0; istep<nsteps; istep++) {

    TIMER_START;

    if(!onlyHeatBath){

      for(j=0; j<overrelaxiter; j++) {

        OverrelaxRed(Ra, Rb, Rc, Ba, Bb, Bc, Jpx, Jpy, Jpz, Jmx, Jmy, Jmz);

        OverrelaxBlue(Ba, Bb, Bc, Ra, Rb, Rc, Jmx, Jmy, Jmz, Jpx, Jpy, Jpz);

      }
    }

    TIMER_STOP;

    ttime+=TIMER_ELAPSED;

    if(!onlyOverrelaxation) {

      //        ComputeRandom(aRand,ZEROFIVE);
      //        ComputeRandom(bRand,ZEROFIVE);
      //        ComputeRandom(cRand,ONEONE);
      //        ComputeRandom(Rand,ZEROONE);

      HeatbathRed(Ra, Rb, Rc, Ba, Bb, Bc, Jpx, Jpy, Jpz, Jmx, Jmy, Jmz);

      HeatbathBlue(Ba, Bb, Bc, Ra, Rb, Rc, Jmx, Jmy, Jmz, Jpx, Jpy, Jpz);

    }


    if(istep && (istep%outfreq)==0) {

      Energy[0]=zero;

      ComputeHEnergy(Ra, Rb, Rc, Ba, Bb, Bc, Jpx, Jpy, Jpz, Jmx, Jmy, Jmz, Energy,TRUE);

      totene=Energy[0];

      printf("Energy per spin at step %d: %f\n",istep,normene*totene);

    }

    if(istep && (istep%dumpfreq)==0) {
      //Dump_Spin(istep);
    }

  }


  printf("Total time for %d iterations of overrelaxation on a %d^3 cube: %f seconds\n", nsteps*overrelaxiter, side, ttime*0.000001);

  printf("Time for spin: %f microseconds\n",ttime/(overrelaxiter*(double)nsteps*V(side)));

  Energy[0]=zero;

  ComputeHEnergy(Ra, Rb, Rc, Ba, Bb, Bc, Jpx, Jpy, Jpz, Jmx, Jmy, Jmz, Energy,TRUE);

  totene=Energy[0];

  printf("Final Energy per spin: %f\n",normene*totene);

  Dump_Spin(nsteps);

  FreeHost();

  return 0;

}


//-------------------------------------------------------------------------------------------------------
// Overrelaxation
//-------------------------------------------------------------------------------------------------------
#define OVERBLOCK factor=two*(as*a[i]+bs*b[i]+cs*c[i])/ \
    (as*as+bs*bs+cs*cs); \
  a[i]=as*factor-a[i]; \
  b[i]=bs*factor-b[i]; \
  c[i]=cs*factor-c[i]; \

#define REDBLOCK as=an[spx]*Jpx[i]+             \
    an[spy]*Jpy[i]+                             \
    an[spz]*Jpz[i]+                             \
    an[smx]*Jmx[i]+                             \
    an[smy]*Jmy[i]+                             \
    an[smz]*Jmz[i];                             \
                                                \
  bs=bn[spx]*Jpx[i]+                            \
    bn[spy]*Jpy[i]+                             \
    bn[spz]*Jpz[i]+                             \
    bn[smx]*Jmx[i]+                             \
    bn[smy]*Jmy[i]+                             \
    bn[smz]*Jmz[i];                             \
                                                \
  cs=cn[spx]*Jpx[i]+                            \
    cn[spy]*Jpy[i]+                             \
    cn[spz]*Jpz[i]+                             \
    cn[smx]*Jmx[i]+                             \
    cn[smy]*Jmy[i]+                             \
    cn[smz]*Jmz[i];                             \
                                                \
  cs+=h;

#define BLUEBLOCK as=an[spx]*Jpx[spx]+                                  \
    an[spy]*Jpy[spy]+                                                   \
    an[spz]*Jpz[spz]+                                                   \
    an[smx]*Jmx[smx]+                                                   \
    an[smy]*Jmy[smy]+                                                   \
    an[smz]*Jmz[smz];                                                   \
                                                                        \
  bs=bn[spx]*Jpx[spx]+                                                  \
    bn[spy]*Jpy[spy]+                                                   \
    bn[spz]*Jpz[spz]+                                                   \
    bn[smx]*Jmx[smx]+                                                   \
    bn[smy]*Jmy[smy]+                                                   \
    bn[smz]*Jmz[smz];                                                   \
                                                                        \
  cs=cn[spx]*Jpx[spx]+                                                  \
    cn[spy]*Jpy[spy]+                                                   \
    cn[spz]*Jpz[spz]+                                                   \
    cn[smx]*Jmx[smx]+                                                   \
    cn[smy]*Jmy[smy]+                                                   \
    cn[smz]*Jmz[smz];                                                   \
                                                                        \
  cs+=h;                                                                \
                                                                        \
                                                                        \
//-------------------------------------------------------------------------------------------------------
// OverrelaxRed
//-------------------------------------------------------------------------------------------------------
void OverrelaxRed(REAL * __restrict__ a, REAL * __restrict__ b, REAL * __restrict__ c,
                  REAL * __restrict__ an, REAL * __restrict__ bn, REAL * __restrict__ cn,
                  REAL * __restrict__ Jpx, REAL * __restrict__ Jpy, REAL * __restrict__ Jpz,
                  REAL * __restrict__ Jmx, REAL * __restrict__ Jmy, REAL * __restrict__ Jmz) {

  int ix, iy, iz;
  REAL as, bs, cs;

  int spx, smx, spy, smy, spz, smz;
  int yoff, zoff, spyoff, smyoff, spzoff, smzoff;

  REAL factor;

  int i, l_maskL=((maskL+1)<<1)-1;
#pragma omp parallel for private(zoff, spzoff, smzoff, yoff, spyoff, smyoff, i, spx, smx, spy, smy, spz, smz, as, bs, cs, factor, ix, iy)
  for(iz=0; iz<side; iz+=2) {

    zoff=(iz<<sdbLsm1);
    spzoff=((iz+1)&l_maskL)<<sdbLsm1;
    smzoff=((iz-1)&l_maskL)<<sdbLsm1;

    for(iy=0; iy<side; iy+=2) { //two lines because of the checkboard!

      yoff=(iy<<sdbLm1);
      spyoff=((iy+1)&l_maskL)<<sdbLm1;
      smyoff=((iy-1)&l_maskL)<<sdbLm1;

      ix=0;
      i=zoff+yoff;
      spx=i;
      smx=i+(Ldb2-1);
      spy=spyoff+zoff;
      smy=smyoff+zoff;
      spz=yoff+spzoff;
      smz=yoff+smzoff;
#if defined(DEBUG)
      printf("iy=%d,iz=%d,i=%d,spx=%d,smx=%d,spy=%d,smy=%d,spz=%d,smz=%d\n",
             iy,iz,i,spx,smx,spy,smy,spz,smz);
#endif
      REDBLOCK;
      OVERBLOCK;

      loopXRed(a, b, c, an, bn, cn,
               Jpx, Jpy, Jpz, Jmx, Jmy, Jmz,
               Ra, Rb, Rc,
               1, Ldb2, yoff+zoff,
               yoff+zoff, yoff+zoff-1, spyoff+zoff, smyoff+zoff, spzoff+yoff, smzoff+yoff);

      yoff=((iy+1)<<sdbLm1);
      spyoff=((iy+2)&l_maskL)<<sdbLm1;
      smyoff=iy<<sdbLm1;
      loopXRed(a, b, c, an, bn, cn,
               Jpx, Jpy, Jpz, Jmx, Jmy, Jmz,
               Ra, Rb, Rc,
               0, Ldb2-1, yoff+zoff,
               yoff+zoff+1, yoff+zoff, spyoff+zoff, smyoff+zoff, spzoff+yoff, smzoff+yoff);

      ix=Ldb2-1;
      i=ix+zoff+yoff;
      spx=zoff+yoff;
      smx=i;
      spy=ix+spyoff+zoff;
      smy=ix+smyoff+zoff;
      spz=ix+yoff+spzoff;
      smz=ix+yoff+smzoff;
#if defined(DEBUG)
      printf("iy=%d,iz=%d,i=%d,spx=%d,smx=%d,spy=%d,smy=%d,spz=%d,smz=%d\n",
             iy+1,iz,i,spx,smx,spy,smy,spz,smz);
#endif
      REDBLOCK;
      OVERBLOCK;

    }
    zoff=((iz+1)<<sdbLsm1);
    spzoff=((iz+2)&l_maskL)<<sdbLsm1;
    smzoff=iz<<sdbLsm1;

    for(iy=0; iy<side; iy+=2) { //two lines because of the checkboard!
      yoff=(iy<<sdbLm1);
      spyoff=((iy+1)&l_maskL)<<sdbLm1;
      smyoff=((iy-1)&l_maskL)<<sdbLm1;

      loopXRed(a, b, c, an, bn, cn,
               Jpx, Jpy, Jpz, Jmx, Jmy, Jmz,
               Ra, Rb, Rc,
               0, Ldb2-1, yoff+zoff,
               yoff+zoff+1, yoff+zoff, spyoff+zoff, smyoff+zoff, spzoff+yoff, smzoff+yoff);

      ix=Ldb2-1;
      i=ix+zoff+yoff;
      spx=zoff+yoff;
      smx=i;
      spy=ix+spyoff+zoff;
      smy=ix+smyoff+zoff;
      spz=ix+yoff+spzoff;
      smz=ix+yoff+smzoff;
#if defined(DEBUG)
      printf("iy=%d,iz=%d,i=%d,spx=%d,smx=%d,spy=%d,smy=%d,spz=%d,smz=%d\n",
             iy,iz+1,i,spx,smx,spy,smy,spz,smz);
#endif
      REDBLOCK;
      OVERBLOCK;

      yoff=((iy+1)<<sdbLm1);
      spyoff=((iy+2)&l_maskL)<<sdbLm1;
      smyoff=iy<<sdbLm1;

      ix=0;
      i=zoff+yoff;
      spx=i;
      smx=i+(Ldb2-1);
      spy=spyoff+zoff;
      smy=smyoff+zoff;
      spz=yoff+spzoff;
      smz=yoff+smzoff;
#if defined(DEBUG)
      printf("iy=%d,iz=%d,i=%d,spx=%d,smx=%d,spy=%d,smy=%d,spz=%d,smz=%d\n",
             iy+1,iz+1,i,spx,smx,spy,smy,spz,smz);
#endif
      REDBLOCK;
      OVERBLOCK;

      loopXRed(a, b, c, an, bn, cn,
               Jpx, Jpy, Jpz, Jmx, Jmy, Jmz,
               Ra, Rb, Rc,
               1, Ldb2, yoff+zoff,
               yoff+zoff, yoff+zoff-1, spyoff+zoff, smyoff+zoff, spzoff+yoff, smzoff+yoff);

    }

  }

}


//-------------------------------------------------------------------------------------------------------
// OverrelaxBlue
//-------------------------------------------------------------------------------------------------------
void OverrelaxBlue(REAL * __restrict__ a, REAL * __restrict__ b, REAL * __restrict__ c,
                   REAL * __restrict__ an, REAL * __restrict__ bn, REAL * __restrict__ cn,
                   REAL * __restrict__ Jpx, REAL * __restrict__ Jpy, REAL * __restrict__ Jpz,
                   REAL * __restrict__ Jmx, REAL * __restrict__ Jmy, REAL * __restrict__ Jmz) {

  int ix, iy, iz;
  REAL as, bs, cs;

  int spx, smx, spy, smy, spz, smz;
  int yoff, zoff, spyoff, smyoff, spzoff, smzoff;

  REAL factor;

  int i, l_maskL=((maskL+1)<<1)-1;
#pragma omp parallel for private(zoff, spzoff, smzoff, yoff, spyoff, smyoff, i, spx, smx, spy, smy, spz, smz, as, bs, cs, factor, ix, iy)
  for(iz=0; iz<side; iz+=2) {

    zoff=(iz<<sdbLsm1);
    spzoff=((iz+1)&l_maskL)<<sdbLsm1;
    smzoff=((iz-1)&l_maskL)<<sdbLsm1;

    for(iy=0; iy<side; iy+=2) { //two lines because of the checkboard!

      yoff=(iy<<sdbLm1);
      spyoff=((iy+1)&l_maskL)<<sdbLm1;
      smyoff=((iy-1)&l_maskL)<<sdbLm1;

      loopXBlue(a, b, c, an, bn, cn,
                Jpx, Jpy, Jpz, Jmx, Jmy, Jmz,
                Ba, Bb, Bc,
                0, Ldb2-1, yoff+zoff,
                yoff+zoff+1, yoff+zoff, spyoff+zoff, smyoff+zoff, spzoff+yoff, smzoff+yoff);

      ix=Ldb2-1;
      i=ix+zoff+yoff;
      spx=zoff+yoff;
      smx=i;
      spy=ix+spyoff+zoff;
      smy=ix+smyoff+zoff;
      spz=ix+yoff+spzoff;
      smz=ix+yoff+smzoff;
#if defined(DEBUG)
      printf("iy=%d,iz=%d,i=%d,spx=%d,smx=%d,spy=%d,smy=%d,spz=%d,smz=%d\n",
             iy+1,iz,i,spx,smx,spy,smy,spz,smz);
#endif
      BLUEBLOCK;
      OVERBLOCK;
      yoff=((iy+1)<<sdbLm1);
      spyoff=((iy+2)&l_maskL)<<sdbLm1;
      smyoff=iy<<sdbLm1;

      ix=0;
      i=zoff+yoff;
      spx=i;
      smx=i+(Ldb2-1);
      spy=spyoff+zoff;
      smy=smyoff+zoff;
      spz=yoff+spzoff;
      smz=yoff+smzoff;
#if defined(DEBUG)
      printf("iy=%d,iz=%d,i=%d,spx=%d,smx=%d,spy=%d,smy=%d,spz=%d,smz=%d\n",
             iy,iz,i,spx,smx,spy,smy,spz,smz);
#endif
      BLUEBLOCK;
      OVERBLOCK;

      loopXBlue(a, b, c, an, bn, cn,
                Jpx, Jpy, Jpz, Jmx, Jmy, Jmz,
                Ba, Bb, Bc,
                1, Ldb2, yoff+zoff,
                yoff+zoff, yoff+zoff-1, spyoff+zoff, smyoff+zoff, spzoff+yoff, smzoff+yoff);

    }

    zoff=((iz+1)<<sdbLsm1);
    spzoff=((iz+2)&l_maskL)<<sdbLsm1;
    smzoff=iz<<sdbLsm1;

    for(iy=0; iy<side; iy+=2) { //two lines because of the checkboard!
      yoff=(iy<<sdbLm1);
      spyoff=((iy+1)&l_maskL)<<sdbLm1;
      smyoff=((iy-1)&l_maskL)<<sdbLm1;

      ix=0;

      i=zoff+yoff;
      spx=i;
      smx=i+(Ldb2-1);
      spy=spyoff+zoff;
      smy=smyoff+zoff;
      spz=yoff+spzoff;
      smz=yoff+smzoff;

#if defined(DEBUG)
      printf("iy=%d,iz=%d,i=%d,spx=%d,smx=%d,spy=%d,smy=%d,spz=%d,smz=%d\n",iy+1,iz+1,i,spx,smx,spy,smy,spz,smz);
#endif

      BLUEBLOCK;
      OVERBLOCK;

      loopXBlue(a, b, c, an, bn, cn,
                Jpx, Jpy, Jpz, Jmx, Jmy, Jmz,
                Ba, Bb, Bc,
                1, Ldb2, yoff+zoff,
                yoff+zoff, yoff+zoff-1, spyoff+zoff, smyoff+zoff, spzoff+yoff, smzoff+yoff);

      yoff=((iy+1)<<sdbLm1);
      spyoff=((iy+2)&l_maskL)<<sdbLm1;
      smyoff=iy<<sdbLm1;

      loopXBlue(a, b, c, an, bn, cn,
                Jpx, Jpy, Jpz, Jmx, Jmy, Jmz,
                Ba, Bb, Bc,
                0, Ldb2-1, yoff+zoff,
                yoff+zoff+1, yoff+zoff, spyoff+zoff, smyoff+zoff, spzoff+yoff, smzoff+yoff);

      ix=Ldb2-1;
      i=ix+zoff+yoff;
      spx=zoff+yoff;
      smx=i;
      spy=ix+spyoff+zoff;
      smy=ix+smyoff+zoff;
      spz=ix+yoff+spzoff;
      smz=ix+yoff+smzoff;
#if defined(DEBUG)
      printf("iy=%d,iz=%d,i=%d,spx=%d,smx=%d,spy=%d,smy=%d,spz=%d,smz=%d\n",
             iy,iz+1,i,spx,smx,spy,smy,spz,smz);
#endif
      BLUEBLOCK;
      OVERBLOCK;

    }
  }

}


//-------------------------------------------------------------------------------------------------------
// ComputeHEnergy
//-------------------------------------------------------------------------------------------------------
void ComputeHEnergy(REAL *a, REAL *b, REAL *c,
                    REAL *an, REAL *bn, REAL *cn,
                    REAL *Jpx, REAL *Jpy, REAL *Jpz,
                    REAL *Jmx, REAL *Jmy, REAL *Jmz,
                    REAL *energy,
                    int even) {

  unsigned int iy, iz;
  int diff;
  int spx, smx, spy, smy, spz, smz;
  REAL e;

  int i;

  for(i=0; i<nspin; i++) {

    iz=i>>sdbLsm1;
    iy=(i - (iz<<sdbLsm1))>>sdbLm1;

    //do the same thing if both $iz and $iy are even or odd
    diff=(((iz&1)^(iy&1))^even);

    spx=(i >> (diff-1)) + (((i-(i&maskL)) + ((i+1)&maskL)) >> (-diff));
    smx=(((i-(i&maskL)) + ((i-1)&maskL)) >> (diff-1)) + (i>>(-diff));
    spy=(i-(i&maskL2))+((i+Ldb2)&maskL2);
    smy=(i-(i&maskL2))+((i-Ldb2)&maskL2);
    spz=(i+Lsdb2) & maskL3;
    smz=(i-Lsdb2) & maskL3;

    e=a[i]*(an[spx]*Jpx[i]+
            an[spy]*Jpy[i]+
            an[spz]*Jpz[i]+
            an[smx]*Jmx[i]+
            an[smy]*Jmy[i]+
            an[smz]*Jmz[i])+
      b[i]*(bn[spx]*Jpx[i]+
            bn[spy]*Jpy[i]+
            bn[spz]*Jpz[i]+
            bn[smx]*Jmx[i]+
            bn[smy]*Jmy[i]+
            bn[smz]*Jmz[i])+
      c[i]*(cn[spx]*Jpx[i]+
            cn[spy]*Jpy[i]+
            cn[spz]*Jpz[i]+
            cn[smx]*Jmx[i]+
            cn[smy]*Jmy[i]+
            cn[smz]*Jmz[i]);

#if defined(DEBUG)
    printf("E=%f\n",e);
#endif
    energy[0]+=e;

  }

}

//-------------------------------------------------------------------------------------------------------
// Heatbath
//-------------------------------------------------------------------------------------------------------
#define HBBLOCK                                                         \
                                                                        \
  anew = RANDZEROFIVE;                                                  \
  bnew = RANDZEROFIVE;                                                  \
                                                                        \
  pl=anew*anew+bnew*bnew;                                               \
                                                                        \
  if(pl<=zerotwentyfive) {                                              \
                                                                        \
    cnew = RANDONEONE;                                                  \
                                                                        \
    hk=SQRT(as*as+bs*bs+cs*cs);                                         \
    bhk=Beta*hk;                                                        \
                                                                        \
    invexpbhk2=EXP(-two*bhk);                                           \
    projection=one+(LOG(invexpbhk2+(one-invexpbhk2)*RANDZEROONE))/bhk;  \
    if(projection<-one) { projection=-one; }                            \
                                                                        \
    invnorm=SQRT((one-cnew*cnew)/pl);                                   \
                                                                        \
    anew*=invnorm;                                                      \
    bnew*=invnorm;                                                      \
                                                                        \
    factor=(as*anew+bs*bnew+cs*cnew)/(hk*hk);                           \
                                                                        \
    anew = anew-factor*as;                                              \
    bnew = bnew-factor*bs;                                              \
    cnew = cnew-factor*cs;                                              \
                                                                        \
    factor=SQRT((one-projection*projection)/(anew*anew+bnew*bnew+cnew*cnew)); \
                                                                        \
    projection/=hk;                                                     \
                                                                        \
    a[i]=as*projection+anew*factor;                                     \
    b[i]=bs*projection+bnew*factor;                                     \
    c[i]=cs*projection+cnew*factor;                                     \
                                                                        \
  }

//#define DEBUG

//-------------------------------------------------------------------------------------------------------
// HeatbathRed
//-------------------------------------------------------------------------------------------------------
void HeatbathRed(REAL * __restrict__ a, REAL * __restrict__ b, REAL * __restrict__ c,
                 REAL * __restrict__ an, REAL * __restrict__ bn, REAL * __restrict__ cn,
                 REAL * __restrict__ Jpx, REAL * __restrict__ Jpy, REAL * __restrict__ Jpz,
                 REAL * __restrict__ Jmx, REAL * __restrict__ Jmy, REAL * __restrict__ Jmz) {

  int ix, iy, iz;
  REAL as, bs, cs;

  REAL anew, bnew, cnew;
  REAL hk, bhk, projection, pl, invnorm;

  REAL invexpbhk2;

  int spx, smx, spy, smy, spz, smz;
  int yoff, zoff, spyoff, smyoff, spzoff, smzoff;

  REAL factor;

  int i, l_maskL=((maskL+1)<<1)-1;

  for(iz=0; iz<side; iz+=2) {

    zoff=(iz<<sdbLsm1);
    spzoff=((iz+1)&l_maskL)<<sdbLsm1;
    smzoff=((iz-1)&l_maskL)<<sdbLsm1;

    for(iy=0; iy<side; iy+=2) { //two lines because of the checkboard!

      yoff=(iy<<sdbLm1);
      spyoff=((iy+1)&l_maskL)<<sdbLm1;
      smyoff=((iy-1)&l_maskL)<<sdbLm1;

      ix=0;
      i=zoff+yoff;
      spx=i;
      smx=i+(Ldb2-1);
      spy=spyoff+zoff;
      smy=smyoff+zoff;
      spz=yoff+spzoff;
      smz=yoff+smzoff;
#if defined(DEBUG)
      printf("iy=%d,iz=%d,i=%d,spx=%d,smx=%d,spy=%d,smy=%d,spz=%d,smz=%d\n",
             iy,iz,i,spx,smx,spy,smy,spz,smz);
#endif
      REDBLOCK;
      HBBLOCK;

      loopXRedHB(a, b, c, an, bn, cn,
                 Jpx, Jpy, Jpz, Jmx, Jmy, Jmz,
                 Ra, Rb, Rc,
                 1, Ldb2, yoff+zoff,
                 yoff+zoff, yoff+zoff-1, spyoff+zoff, smyoff+zoff, spzoff+yoff, smzoff+yoff);

      yoff=((iy+1)<<sdbLm1);
      spyoff=((iy+2)&l_maskL)<<sdbLm1;
      smyoff=iy<<sdbLm1;
      loopXRedHB(a, b, c, an, bn, cn,
                 Jpx, Jpy, Jpz, Jmx, Jmy, Jmz,
                 Ra, Rb, Rc,
                 0, Ldb2-1, yoff+zoff,
                 yoff+zoff+1, yoff+zoff, spyoff+zoff, smyoff+zoff, spzoff+yoff, smzoff+yoff);

      ix=Ldb2-1;
      i=ix+zoff+yoff;
      spx=zoff+yoff;
      smx=i;
      spy=ix+spyoff+zoff;
      smy=ix+smyoff+zoff;
      spz=ix+yoff+spzoff;
      smz=ix+yoff+smzoff;
#if defined(DEBUG)
      printf("iy=%d,iz=%d,i=%d,spx=%d,smx=%d,spy=%d,smy=%d,spz=%d,smz=%d\n",
             iy+1,iz,i,spx,smx,spy,smy,spz,smz);
#endif
      REDBLOCK;
      HBBLOCK;

    }
    zoff=((iz+1)<<sdbLsm1);
    spzoff=((iz+2)&l_maskL)<<sdbLsm1;
    smzoff=iz<<sdbLsm1;

    for(iy=0; iy<side; iy+=2) { //two lines because of the checkboard!
      yoff=(iy<<sdbLm1);
      spyoff=((iy+1)&l_maskL)<<sdbLm1;
      smyoff=((iy-1)&l_maskL)<<sdbLm1;

      loopXRedHB(a, b, c, an, bn, cn,
                 Jpx, Jpy, Jpz, Jmx, Jmy, Jmz,
                 Ra, Rb, Rc,
                 0, Ldb2-1, yoff+zoff,
                 yoff+zoff+1, yoff+zoff, spyoff+zoff, smyoff+zoff, spzoff+yoff, smzoff+yoff);

      ix=Ldb2-1;
      i=ix+zoff+yoff;
      spx=zoff+yoff;
      smx=i;
      spy=ix+spyoff+zoff;
      smy=ix+smyoff+zoff;
      spz=ix+yoff+spzoff;
      smz=ix+yoff+smzoff;
#if defined(DEBUG)
      printf("iy=%d,iz=%d,i=%d,spx=%d,smx=%d,spy=%d,smy=%d,spz=%d,smz=%d\n",
             iy,iz+1,i,spx,smx,spy,smy,spz,smz);
#endif
      REDBLOCK;
      HBBLOCK;

      yoff=((iy+1)<<sdbLm1);
      spyoff=((iy+2)&l_maskL)<<sdbLm1;
      smyoff=iy<<sdbLm1;

      ix=0;
      i=zoff+yoff;
      spx=i;
      smx=i+(Ldb2-1);
      spy=spyoff+zoff;
      smy=smyoff+zoff;
      spz=yoff+spzoff;
      smz=yoff+smzoff;
#if defined(DEBUG)
      printf("iy=%d,iz=%d,i=%d,spx=%d,smx=%d,spy=%d,smy=%d,spz=%d,smz=%d\n",
             iy+1,iz+1,i,spx,smx,spy,smy,spz,smz);
#endif
      REDBLOCK;
      HBBLOCK;

      loopXRedHB(a, b, c, an, bn, cn,
                 Jpx, Jpy, Jpz, Jmx, Jmy, Jmz,
                 Ra, Rb, Rc,
                 1, Ldb2, yoff+zoff,
                 yoff+zoff, yoff+zoff-1, spyoff+zoff, smyoff+zoff, spzoff+yoff, smzoff+yoff);

    }

  }

}

//-------------------------------------------------------------------------------------------------------
// HeatbathBlue
//-------------------------------------------------------------------------------------------------------
void HeatbathBlue(REAL * __restrict__ a, REAL * __restrict__ b, REAL * __restrict__ c,
                  REAL * __restrict__ an, REAL * __restrict__ bn, REAL * __restrict__ cn,
                  REAL * __restrict__ Jpx, REAL * __restrict__ Jpy, REAL * __restrict__ Jpz,
                  REAL * __restrict__ Jmx, REAL * __restrict__ Jmy, REAL * __restrict__ Jmz) {

  int ix, iy, iz;
  REAL as, bs, cs;

  REAL anew, bnew, cnew;
  REAL hk, bhk, projection, pl, invnorm;

  REAL invexpbhk2;

  int spx, smx, spy, smy, spz, smz;
  int yoff, zoff, spyoff, smyoff, spzoff, smzoff;

  REAL factor;

  int i, l_maskL=((maskL+1)<<1)-1;

  for(iz=0; iz<side; iz+=2) {


    zoff=(iz<<sdbLsm1);
    spzoff=((iz+1)&l_maskL)<<sdbLsm1;
    smzoff=((iz-1)&l_maskL)<<sdbLsm1;

    for(iy=0; iy<side; iy+=2) { //two lines because of the checkboard!


      yoff=(iy<<sdbLm1);
      spyoff=((iy+1)&l_maskL)<<sdbLm1;
      smyoff=((iy-1)&l_maskL)<<sdbLm1;

      loopXBlueHB(a, b, c, an, bn, cn,
                  Jpx, Jpy, Jpz, Jmx, Jmy, Jmz,
                  Ba, Bb, Bc,
                  0, Ldb2-1, yoff+zoff,
                  yoff+zoff+1, yoff+zoff, spyoff+zoff, smyoff+zoff, spzoff+yoff, smzoff+yoff);

      ix=Ldb2-1;
      i=ix+zoff+yoff;
      spx=zoff+yoff;
      smx=i;
      spy=ix+spyoff+zoff;
      smy=ix+smyoff+zoff;
      spz=ix+yoff+spzoff;
      smz=ix+yoff+smzoff;
#if defined(DEBUG)
      printf("iy=%d,iz=%d,i=%d,spx=%d,smx=%d,spy=%d,smy=%d,spz=%d,smz=%d\n",
             iy+1,iz,i,spx,smx,spy,smy,spz,smz);
#endif
      BLUEBLOCK;
      HBBLOCK;
      yoff=((iy+1)<<sdbLm1);
      spyoff=((iy+2)&l_maskL)<<sdbLm1;
      smyoff=iy<<sdbLm1;

      ix=0;
      i=zoff+yoff;
      spx=i;
      smx=i+(Ldb2-1);
      spy=spyoff+zoff;
      smy=smyoff+zoff;
      spz=yoff+spzoff;
      smz=yoff+smzoff;
#if defined(DEBUG)
      printf("iy=%d,iz=%d,i=%d,spx=%d,smx=%d,spy=%d,smy=%d,spz=%d,smz=%d\n",
             iy,iz,i,spx,smx,spy,smy,spz,smz);
#endif
      BLUEBLOCK;
      HBBLOCK;

      loopXBlueHB(a, b, c, an, bn, cn,
                  Jpx, Jpy, Jpz, Jmx, Jmy, Jmz,
                  Ba, Bb, Bc,
                  1, Ldb2, yoff+zoff,
                  yoff+zoff, yoff+zoff-1, spyoff+zoff, smyoff+zoff, spzoff+yoff, smzoff+yoff);

    }

    zoff=((iz+1)<<sdbLsm1);
    spzoff=((iz+2)&l_maskL)<<sdbLsm1;
    smzoff=iz<<sdbLsm1;

    for(iy=0; iy<side; iy+=2) { //two lines because of the checkboard!
      yoff=(iy<<sdbLm1);
      spyoff=((iy+1)&l_maskL)<<sdbLm1;
      smyoff=((iy-1)&l_maskL)<<sdbLm1;

      ix=0;

      i=zoff+yoff;
      spx=i;
      smx=i+(Ldb2-1);
      spy=spyoff+zoff;
      smy=smyoff+zoff;
      spz=yoff+spzoff;
      smz=yoff+smzoff;

#if defined(DEBUG)
      printf("iy=%d,iz=%d,i=%d,spx=%d,smx=%d,spy=%d,smy=%d,spz=%d,smz=%d\n",iy+1,iz+1,i,spx,smx,spy,smy,spz,smz);
#endif

      BLUEBLOCK;
      HBBLOCK;

      loopXBlueHB(a, b, c, an, bn, cn,
                  Jpx, Jpy, Jpz, Jmx, Jmy, Jmz,
                  Ba, Bb, Bc,
                  1, Ldb2, yoff+zoff,
                  yoff+zoff, yoff+zoff-1, spyoff+zoff, smyoff+zoff, spzoff+yoff, smzoff+yoff);

      yoff=((iy+1)<<sdbLm1);
      spyoff=((iy+2)&l_maskL)<<sdbLm1;
      smyoff=iy<<sdbLm1;

      loopXBlueHB(a, b, c, an, bn, cn,
                  Jpx, Jpy, Jpz, Jmx, Jmy, Jmz,
                  Ba, Bb, Bc,
                  0, Ldb2-1, yoff+zoff,
                  yoff+zoff+1, yoff+zoff, spyoff+zoff, smyoff+zoff, spzoff+yoff, smzoff+yoff);

      ix=Ldb2-1;
      i=ix+zoff+yoff;
      spx=zoff+yoff;
      smx=i;
      spy=ix+spyoff+zoff;
      smy=ix+smyoff+zoff;
      spz=ix+yoff+spzoff;
      smz=ix+yoff+smzoff;
#if defined(DEBUG)
      printf("iy=%d,iz=%d,i=%d,spx=%d,smx=%d,spy=%d,smy=%d,spz=%d,smz=%d\n",
             iy,iz+1,i,spx,smx,spy,smy,spz,smz);
#endif
      BLUEBLOCK;
      HBBLOCK;

    }
  }

}
